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THESIS DISCLAIMER 



The reader is cautioned that computer programs developed in 
this research may not have been exercised for all cases of 
interest. While every effort has been made within the time 
available to ensure that the programs are free of computational and 
logic errors, they cannot be considered validated. Any application 
of these programs without additional verification is at the risk 
of the user. 
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I. 



INTRODUCTION 



In many radar and remote sensing operations, successive 
image frames may contain a moving object of interest. 
Identification of this object may be complicated by image 
degradation resulting from the presence of background noise. 
The idea of using an extended Kalman filter (EKF) to enhance 
the quality of the image frame and provide an estimate of the 
object velocity was proposed by Burl [Ref. 1] . 

In this thesis, a moving image model is defined as a 
series of image frames containing a moving object. An image 
frame consists of a two-dimensional array where the individual 
array elements are defined as pixels. Each pixel is given a 
numerical value based on the intensity of the image at that 
point. These values are either stored in a computer or 
projected on a video display terminal. 

The EKF requires a model that describes the evolution of 
the image frames in time. This spatial domain model is 
derived in Chapter II. A shift operator in the model 
describes the motion of the object as it traverses the image 
frame and is a function of the object's velocity. A 
prohibitive number of calculations would be necessary if the 
EKF were implemented in the spatial domain. [Ref. 1] 

Chapter II further describes the moving image model after 
it has been transformed to the spatial frequency domain. The 
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two-dimensional discrete Fourier transform (DFT) transforms 
these two-dimensional image frames from the spatial domain to 
the spatial frequency domain. By transforming the image frame 
to the spatial frequency domain, the EKF can analyze each of 
the spatial frequencies separately. The moving object is 
evidenced by a phase shift in each of the spatial frequencies. 
The implementation of the EKF in the spatial frequency domain 
reduces the number of calculations. The number of 
calculations can be further reduced by limiting the number of 
spatial frequencies that are analyzed. [Ref. 1] 

Chapter III presents the extended Kalman filter equations. 
As the velocity estimate error approaches zero, the EKF is 
shown to converge to a parallel set of third-order extended 
Kalman filters. This parallel structure is referred to as the 
modified extended Kalman filter (MEKF) . [Ref. 1] 

Chapter IV discusses the results of the computer 
simulations that were run to record the performance of the 
EKF under a variety of noise conditions and background 
environments. Conclusions are presented in Chapter V. 
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II. THE MOVING IMAGE MODEL 



A. THE MOVING IMAGE MODEL IN THE SPATIAL DOMAIN 

Successive image frames containing a moving object can be 
modeled by a nonlinear state equation in the spatial domain 
as, 



x(k+l) 



f(k+l) 




S(v) 0 




f(k) 




Cj(k) 


v(k+l) 




0 I 




v(k) 


1 





= a(v)x(k) + C(k) 



( 2 . 1 ) 



An image frame has the dimensions of Ni by N 2 pixels. The 
amplitude of each pixel, when stacked into a column, forms 
the vector f, where f € The velocity vector, v, 

consists of a two-element column vector describing the 
horizontal and vertical velocities of the moving object in 
terms of the number of pixels per iteration. The vectors f 
and v combine to form the state vector x, where x e 
S(v) is a two-dimensional shift operator with the magnitude 
and direction being a function of the velocity vector, v. 
The plant noise vector, f, is composed of the image plant 
noise, and the velocity plant noise, The plant noise 

is assumed to be a zero-mean, white, Gaussian random process 
with a covariance matrix defined as 
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( 2 . 2 ) 



Qn = E 



[ 




t 5:, 




a^fl 0 
0 o\l 




where E[*] is the expectation operator and I is the identity 
operator. 

The measured image, y(k) , is defined by the corresponding 
measurement equation, 

y(k) = f(k) + v(k) 

= [ I 0 ] x(k) + v(k) 

= c x(k) + v(k) , (2.3) 

, (ninp) , (ninp) , , . 

where ye/? ' and u g R ‘ ^ . The zero mean, white, Gaussian 
random process, i/(k), has a known covariance matrix 

R,, = E [ v(k) v^(k) ] = o\l . (2.4) 



B. DERIVATION OF THE SHIFT OPERATOR IN THE SPATIAL FREQUENCY 
DOMAIN 

A image frame is a two-dimensional, finite-extent sequence 
defined as f(ni,n 2 ) where 0 < n^ < Ni-1 and 0 < n 2 < N 2 -I. The 
sequence f(ni,n 2 ) is transformed from the spatial domain to the 
spatial frequency domain using the two-dimensional discrete 
Fourier transform (DFT) , 



F(k,,k2) 



N|-l N2-I 

I I 

nj=0 1)2 = 0 



f (n,,n2)exp 





(2.5) 
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where 0 < ki < Ni-1 and 0 < k 2 < N 2 -I. 

The shift operator in the spatial frequency domain is 
derived from the circular shift property of the discrete 
Fourier transform [Ref. 2:p. 67]. The circular shift property 
is defined as 

x( ((n,-m,))^^, ((n^-m^))^^ ) x(k,,k 2 ) 

♦-> D(m, ,m 2 )X(k, ,k 2 ) , (2.6) 

-J— 

where is the modulus operator, = e and D(m[,m 2 ) 

is the transformed shift operator. Equation 2.6 indicates 
that a circular shift in the spatial domain results in a phase 
shift in the spatial frequency domain. 

The original image frame, x(ni,n 2 ), is assumed to be a 
periodic two-dimensional sequence with the fundamental period 
equal to the frame dimensions. The object moving across an 
image frame then describes a circular shift. [Ref. 2:pp. 61- 
62]. This implies that, as the object moves off one edge of 
the screen, it reenters the screen from the opposite side. 
This assumption was necessary because the processing 
capabilities of the computer used to execute the EKF algorithm 
limited the image frame dimensions. Image frames of finite 
extent, such as a radar screen or video display, do not 
require this property because the moving object's size will 
likely be disproportionally smaller than the screen 
dimensions. 
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Derivation of D(mi,m2) begins by introducing a circular 
shift into the original image frame, f(ni,n2). Equation 2.5 
now becomes , 



F'(ki,k 2 ) 



N,-l N2-I 

i S f (n,-mi ,n2-m2) 

Oj H 2 = 0 



2 irnj k j 



2 ^ 2 ^-^ 

N2 J 



( 2 . 7 ) 



where F'(ki,k2) is the circular-shifted, transformed image 
frame. The amount of pixel movement is described by mi for 
the horizontal shift and m2 for the vertical shift. 

Letting li = ni-mi and I2 = n2~m2. Equation 2.7 becomes 



F'(k,,k2) = S 



N,-mi-l N 2 -m 2 -l 



l,r -m, I2 = -n>2 



f(l,,l2) e 



^ 2it( 1 j )k j 2ir( I 2 + 0^2 )k2 

I " ^2 J 



. (2.8) 



Rearranging terms. Equation 2.8 reduces to 

^ 2 Tim I k I 2 ^^2 ^2 

F'(k,,k2) = e ^ ^ F(k,,k2) 

n) 1 k 1 nin k 

= Wk, ‘ 2 2 F(k,,k2) 

= D(mi,m2) F(k],k2) , (2.9) 

where F(ki,k2) is the transformed image frame without the 
circular shift. This proves that, when an object moves across 
an image frame, the motion is described by a phase shift in 
the spatial frequency domain, thereby agreeing with the 
circular shift property (Eqn. 2 . 6 ). 
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The notation is simplified by defining the frequency 



vector, w, as 





” 27tk, ” 


r 1 - 


N. 


L J 


27tk2 




L N2 J 



and the velocity vector as 



V 






( 2 . 10 ) 



( 2 . 11 ) 



Equation 2 . 9 now becomes 

F' (k,,k2) = F(k,,k2) 

= D(m,,m2) F(k,,k2) 

= D(v) F(k,,k 2 ) , (2.12) 

where T is the transpose operation. 

The example contained in Appendix A further demonstrates 
the use of the phase shift operator in reconstructing image 
frames in the spatial frequency domain. 



C. THE MOVING IMAGE MODEL IN THE SPATIAL FREQUENCY DOMAIN 

The spatial domain model described by Equations 2 . 1 and 

2.3 can now be transformed to the spatial frequency domain. 

The transformation is accomplished by using the two- 

dimensional discrete Fourier transform, 

N|-i N 2-1 p fnk nk'^n 

F(k,,k2) = I S f(n,,n2)exp -j27t ,(2.13) 

ri|=0 n2=0 l_ v I 2 y J 

where 0 < ki < Ni~l and 0 < k 2 < N 2 -I. 
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The vector, f(k), and its corresponding plant noise, f 



are transformed to the spatial frequency domain. 





- Re(F, (k)) 




Im(F,(k)) 


F(k) = F{f(k) ) = 


Re I 






Im| 








v(k) 



f / 



(2.14) 



Z(k) = } , (2.15) 
where f{*} is the two-dimensional discrete Fourier transform. 
(It should be noted that the index k refers to the sample at 
time k and is not related to the frequency vector.) Only half 
of the spatial frequencies need to stored in F(k) due to the 
conjugate symmetry that results from the transformation of a 
real image [Ref. 3;p. 45]. 

The state equation is defined in the spatial frequency 
domain as. 



X(k+1) = A(v)X(k) + Z(k) 



rF(k+l)-| ^ FD(v) 0-irF(k)-| ^ rzF(^)"| 

Lv(k+i)J L 0 iJLv(k)J J ' 



(2.16) 



where X e rj-^e transformed shift operator, D(v) , is a 

diagonal matrix where the diagonal terms, D,- (v) , are the 
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transformed shift operators for a specific spatial frequency. 
Because the real and imaginary parts of the spatial frequency 
vector are stacked into a column vector, D,- (v) is expressed 
in terms of its trigonometric relations, 

F,(k+1) = D,(v)Fj(k) 



rRe(F, 


(k+1)) 1 


r cos(w,V) 


sin(WjV) ”1 


rRe(F, 


(k)) 1 


L 


(k+D) J" 


-sin(WjV) 


COS(WjV) J 


L IinCFi 


(k)) J 



where i 



0 , 1 ,.. 



(Ni-1) (N^-l) 
2 



The covariance matrix for the image plant noise, Z, is. 



= E [z Z^] 




0 

0^ I 




(2.18) 



where E[*] is the expectation operator and I is the identity 
operator. Because the discrete Fourier transform is a linear 
operation on Z remains a zero-mean, Gaussian random 
process . 

The corresponding measurement equation in the spatial 
frequency domain becomes. 



Y(k) = F(k) + N(k) 

= [ I 0 ] X(k) + N(k) 

= C X(k) + N(k) (2.19) 
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where Y(k) = f {y(k) ) and N(k) = F{i/(k)). The state equations 
given by Equations 2.16 and 2.19 are the basis for 
implementing the EKF in the spatial frequency domain. 
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III. THE EXTENDED KAIiMAN FILTER 



A. GENERAL 

In this chapter, the extended Kalman filter (EKF) is 
obtained from the state equations that were derived for the 
moving image model. Equations 2.16 and 2.19 are used to 
implement the EKF in the spatial frequency domain. The 
modified extended Kalman filter (MEKF) is based on the 
diagonal properties of the transformed shift operator, D(v) , 
and covariance matrices. The parallel structure of the MEKF 
is formed in the limit as the velocity estimate approaches 
the actual velocity. 

B. THE EXTENDED KALMAN FILTER 

The extended Kalman filter (EKF) algorithm is used when 
the state equations are linearized about each new estimate as 
they become available [Ref. 4: p.l58 ]. The EKF uses the 

spatial frequency domain state equations, Eqns 2.16 and 2.19, 
for the moving image model. 

The linearized state equation is given as 

6X(k+l) = A(XJ6X(k) + Z(k) , (3.1) 

where 5 X (k) =X (k) -X q. Here, Xq is shorthand notation for the 
current estimate X(k|k) [Ref. 4: p.l58] and A(Xq) is the 

Jacobian of the nonlinear state equations about Xq. This 
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means that the z 7 component of A(Xo) is the partial derivative 
with respect to the 7 th state of the zth component of the state 
equation (Eqn. 2.16) and A(Xq) is then composed as follows: 



A(XJ = diagonal [A, (X,) ] 



i = 1,2, .. . (3.2) 



N.N, 



cos(a) sin(a) [-f , (k) sin(a) +f 2 (k) cos (a) ” 

-sin(a) cos(a) [-f, (k) cos (a) -f 2 (k) sin (a) /(3.3) 



0 



0 



1 



where a = w^v, f, (k) = Re[F,(k)] and f 2 (k) = Im[Fj(k)]. A 
further discussion of the linearization techniques used by 
EKFs is contained in Ref. 4, p.l54, and Ref. 5, p.l94. 

The EKF equations for the moving image model are 
summarized below. 

Prediction 

State prediction 



X(k+l|k) = A(v) X(k|k) 



(3.4) 



Covariance prediction 



P(k+l|k) = A(XJX(k|k)P(k|k)A^XJX(k|k) + Q 22 (3.5) 



Innovation 



e(k+l) = Y(k+1) Y(k+l|k) 

= Y(k+1) - C X(k+l|k) 



(3.6) 



Kalman gain 



G(k+1) = P(k+l|k)C^[C P(k+l|k) d + Rnn]‘‘ (3-V) 
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Correction 



State correction 

X(k+l|k+l) = X(k+l|k)+G(k+l) [Y(k+l)-Y(k+l|k) ] (3.8) 
Covariance correction 

P(k+l|k+l) = [ I - G(k+l)C]P(k+l|k) (3.9) 

The notation k+l|k reads as "at time k+1, given data up 
through time k.” 

B. MODELING PROBLEMS 

When the Kalman filter is implemented, the operation of 
the filter can be degraded by modeling errors resulting from 
linearization around the imperfectly estimated state. One 
possible method of reducing this divergence is to add 
fictitious plant noise to the system model. [Ref. 6: p.280] 
The system dynamics matrix, A(Vo+A) , is assumed to have 
model errors which is represented by A in Eqn. 3.10 and is 
expressed as a first-order Taylor series expansion. The 
first-order term is then combined with the actual plant noise 
vector, w(k) , to form the plant noise vector, f. 



X(k+1) = A(v^+A)X(k) + w(k) 



= A(v)X(k) + 



[ 



3A(v) 

3v 



v=v 



0 



A X(k) + w 




= A(v)X(k) + C(k) 



(3.10) 



where A denotes the model error and w(k) is the actual plant 
noise vector. Including the model error, as part of the plant 
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noise vector, prevents new measurements from being weighted 
too lightly because the Kalman gains tend to zero as t — * oo . 
This ensures that the model continues to track the system as 
new measurements are taken. 

C. THE MODIFIED EXTENDED KAIMAN FILTER 

As mentioned earlier, the diagonal properties of the 
transformed shift operator, D(v) , and the covariance matrices 
are the basis for the parallel structure shown in Figure 3.1. 

This parallel structure is referred to as the modified 
extended Kalman filter. Reference 1 outlines how the EKF 
converges to the MEKF as the velocity estimate approaches the 
actual velocity. 

Practical implementation of the MEKF suggests that each 
filter in the parallel bank be dedicated to a specific spatial 
frequency. Each individual filter is referred to as a single 
frequency extended Kalman filter (SFEKF) . 

The state vector for a specific spatial frequency is 
defined as: 





“X, ,(k) “ 




-Re(F,(k)) - 


X,(k) = 


X,.2(k) 


= 


Im(F,(k)) 




_ X, 3(k) _ 




(i>iV 



where the first two states, X,- •|(k) and X,- 2 (h) , are the Fourier 
coefficients associated with a specific spatial frequency and 
the third state, X,- 3 (k), is the dot product of the spatial 
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frequency vector and the velocity vector. The reason for 
making X,- 3 (k) a dot product is because the velocity vector is 
not directly observable for a specific spatial frequency. 
Even though the Kalman filter can generate estimates without 
the state being observable, estimation of the unobservable 
portion of the state relies on a priori information and not 
on the measurements [Ref. 1]. Simplification of the 
individual filters results because only Xj^ 3 (k) is estimated. 
The a priori information is later combined with all the 
estimates of X,-^ 3 (k) to yield a final velocity estimate. The 
final velocity estimate is calculated using the weighted least 
squares algorithm given by: 



v(k) = (w^z‘'o))''oj^z''x 3 ,| 3 (k) , (3.12) 

where w was defined by Eqn. 2.10 and 



1 . 3.3 






z ‘ = 



N1N2 



.3,3 






(3.13) 



Pj 3 3 is the 3,3 component of the /th estimation error 
covariance matrix computed by each of the SFEKF's. 
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Figure 3.1 The Modified Extended Kalman Filter 



16 



The state equation for the SFEKF is 



Xi(k+1) = A,(v)X,(k) + Z,(k) 
rD,(v) 0 -| 

= L ‘o 1 J 2,(k) 

and Zj has the known covariance matrix 
= ECZ,Z/] - 



0 0 
0 0 
0 0 llw, 11^0^2 



(3.14) 



(3.15) 



The measurement equation and corresponding measurement 
noise covariance matrix for the SFEKF are similarly derived 
for each spatial frequency and have the form shown earlier in 
Eqn. 2.19. 

The linearized state equation is given by 

8X,(k+l) = A,(X^^)8X,(k) + Z,(k) , (3.16) 

where Aj(Xg^)is the Jacobian for a specific spatial frequency 
and is linearized about X . 

The extended Kalman filter equations, Eqns. 3.3 through 
3.8, are similarly modified to accommodate the SFEKF model. 

D. THE REDUCED ORDER EXTENDED KALMAN FILTER 

Data compression of images is achieved by truncating 
spatial frequencies [Ref. 1] . The parallel structure of the 
MEKF can easily accommodate this truncation and suggests a 
reduced-order filter. Even though this will lead to a 
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deterioration of the output image quality, the degrading 
effects of the image noise are simultaneously reduced. 

The actual selection of,, which spatial frequencies to 
truncate will depend upon the distribution of the energy in 
the signal and noise level. A priori knowledge of the spatial 
frequency content would assist in the selection of which 
spatial frequencies are to be analyzed when filtering a series 
of image frames. 
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IV. SIMULATIONS 



The MATLAB program written for this thesis simulated an 
image frame with the dimensions of 32x32 pixels and is 
contained in Appendix C. The two types of moving objects 
employed for test purposes were an 8x8 pixel square of unity 
amplitude (Figure. 4.1) and an 8x8 pixel pyramid with 
amplitude intervals of 0.25 (Figure 4.2). Additionally, tests 
were conducted using a fixed a checkerboard background. 
Figures 4.3 and 4.4 depict the test objects moving across a 
checkerboard pattern. The individual checkers were 8x8 pixel 
squares with an amplitude of 0.25. 

A. PERIODICITY OF THE TEST OBJECTS 

An interesting challenge had to be overcome because of 
the small dimensions of the image frame. The circular shift 
property of the discrete Fourier transform allowed the 
computer simulation to be carried out over several hundred 
iterations. The original 32x32 pixel image frame used in the 
simulations was assumed to be a periodic two-dimensional 
sequence. The moving object quickly traversed the image frame 
because of the small dimensions. By assuming periodicity, 
test objects that went off one edge of the frame would 
reappear on the opposite edge. This assumption allowed 
simulations to be run over several hundred iterations because 
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Figure 4 . 1 
amplitude. 



An 8x8 pixel, square object of unity 




Figure 4.2 An 8x8 pixel, pyramid object with amplitude 
intervals of 0.25. 
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Figure 4.3 An 8x8 pixel, square object of unity 
amplitude moving across a checkerboard background where 
each of the individual checkers has an amplitude of 
0.25. 




Figure 4.4 An 8x8 pixel, pyramid object with amplitude 
intervals of 0.25 moving across a checkerboard 
background where the individual checkers have an 
amplitude of 0.25. 
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it appeared to the extended Kalman filter that the image frame 
was of a larger dimension. 

B. DISCRETE FOURIER TRANSFORM OF MOVING IMAGES 

A description of how the transformed image frame appears 
in the spatial frequency domain will help demonstrate why the 
extended Kalman filter was implemented in the spatial 
frequency domain. 

As an example, the square object depicted in Figure 4.1 
is transformed to the spatial frequency domain and the 
resulting magnitude plot is shown in Figure 4.5. Just as a 
one-dimensional square wave transforms to a one-dimensional 
sine function, the two-dimensional square object transforms 
to a two-dimensional sine function. 

What is not obvious is that as the square object is moved 
across the image frame, a linear phase shift results in each 
of the spatial frequencies; the magnitude plot remains the 
same. The amount of phase shift is frequency dependent as 
was shown in Equation 2.21. 

Figure 4.6 shows the linear phase shift that occurs in 
the spatial frequencies F(1,0) and F(0,1). When the square 
object is moved horizontally 1 pixel per iteration, a phase 
shift of -11.25° occurs in F(1,0). Likewise, when the square 
image is simultaneously shifted vertically 2 pixels per 
iteration, F(0,1) experiences a 22.5° phase shift per 
iteration. 
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Figure 4.5 Magnitude plot of the transformed image 
frame depicted in Fig. 4.1. 
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Figure 4.6 A linear phase shift of -11.25° per iteration 
occurs in F(1,0) (dashed line) and 22.5° per iteration for 
F(0,1) (solid line) when the square object is moved 1 
pixel horizontally and 2 pixels vertically per iteration. 
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Appendix B further discusses the idea of how objects 
moving at a constant velocity in a no-noise environment create 
a linear phase shift in the spatial frequency domain. Also 
shown is the effect of adding white noise to the image frame 
before it is transformed to the spatial frequency domain. 

C. LOWPASS FILTERING OF THE TRANSFORMED IMAGES 

After the test objects were transformed to the spatial 
frequency domain, the specific spatial frequencies were 
analyzed. Most of the energy for the two test objects was 
centered about the dc component which was located at the 
center of the transformed array. 

Because most of the energy was found near the dc 
component, the high frequency components were truncated or 
filtered out using a lowpass filter, 

Y(ki,k2) = X(k,,k2)H(ki,k2) (4.1) 

where Y(ki,k 2 ) is the filtered image, X(ki,k 2 ) is the 
transformed image, and H(ki,k 2 ) is the lowpass filter. 

As an example. Figure 4.1 was transformed to the spatial 
frequency domain where a lowpass filter was applied. This 
operation truncated the high frequency components. After 
transformation back to the spatial domain. Figure 4.7 shows 
the resulting filtered image frame. This filtered image frame 
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after being filtered in the spatial frequency domain. 
The sharp edges are lost due to the loss of the high 
frequency components. 
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will be used in finding the error to calculate the Frobenius 
norm, which is discussed in the next section. 



D. PERFORMANCE MEASURES OF THE EXTENDED KALMAN FILTER 

After the specified number of iterations had been 
completed for each simulation, the final, filtered image frame 
was transformed back to the spatial domain. This filtered 
image frame could then be visually compared with the no-noise 
image frame to see how effective the extended Kalman filter 
performed. 

An additional performance measure involved analytically 
calculating the Frobenius norm throughout the simulation. 
The Frobenius norm is a scalar which is computed by squaring 
each of the terms in a matrix, summing each of the squares, 
and then taking the square root of the final sum (Eqn. 4.1) 
[Ref. 7]. 



Two Frobenius norms were calculated for each of the 
simulations. The first Frobenius norm was calculated from 
the error matrix that was found by subtracting the current 
estimate of the image frame from the current, no-noise image 
frame. Because the extended Kalman filter only analyzed 
specific spatial frequencies, the current no-noise image frame 




(4.1) 
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was lowpass filtered so that only those spatial frequencies 
analyzed by the extended Kalman filter remained. As was shown 
earlier, Figure 4.1 is an example of a current, no-noise image 
while its lowpass filtered counterpart is depicted in Figure 
4.7. The current estimate of the image frame was then 
subtracted from the lowpass filtered image frame and the 
second Frobenius norm was calculated for this error matrix. 

E. SIMULATION RESULTS 

The first set of computer simulations had a square test 
object traversing the image frame. The extended Kalman filter 
was evaluated under a variety of test conditions which are 
summarized in Table 4.1. 

The two types of backgrounds were a homogeneous background 
with zero intensity and a checkerboard background. The 
velocity in Table 4.1 is described in parentheses where the 
first number is the horizontal velocity and the second number 
is the vertical velocity. Finally, the standard deviation of 
the white noise describes the noise level that was introduced 
to the image frame during the successive iterations. 
Additionally, to ensure that the same random noise sequence 
was used for each of the computer simulations, a random number 
seed was used to initilize the random number generator. 

Figures 4.8 through 4.47 are the computer simulation 
results using the square test object. The no-noise image 
frame, that the extended Kalman filter is estimating, is 
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frame that the extended Kalman filter is estimating is 
presented first. The extended Kalman filter estimate is then 
depicted; followed by the horizontal and vertical velocity 
estimates and the Frobenius norm plots. 

TABLE 4.1 SIMULATION TEST CONDITIONS (SQUARE OBJECT) 



Simulation 


Background 


Velocity 


Noise level 


1 


H 


(1.0, 2.0) 


1.0 


2 


C 


(1. 0, 2. 0) 


1.0 


3 


H 


(1.0, 2.0) 


2.0 


4 


C 


(1.0, 2.0) 


2.0 


5 


H 


(1.0, 2.0) 


4.0 


6 


C 


(1.0, 2.0) 


4.0 


7 


H 


(2. 5, 1.5) 


1.0 


8 


C 


(2.5, 1.5) 


1.0 


9 


H 


(1.0, 0.5) 


1.0 


10 


C 


H 

• 

O 

O 

• 


1.0 



H - homogeneous background, C - checkerboard background 



The results of the first two simulations are shown in 
Figures 4.8 through 4.15. The noise level for both 
simulations had a standard deviation of 1.0. The object was 
moved horizontally l pixel per iteration and simultaneously 
moved vertically 2 pixels per iteration. Even with the 
checkerboard background, the EKF was able to rapidly estimate 
the object velocity for both simulations. The Frobenius norm 
correspondingly reached steady state as the velocity estimates 
approached the actual object velocity. 
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The next two simulations are depicted in Figures 4 . 16 
through 4.23. The standard deviation of the noise was 
increased to 2.0. Figures -^4.16 and 4.20 show the square 
object position after 95 iterations. These two figures 
demonstrate the assumption that the image frame used in these 
simulations was a two-dimensional, periodic sequence. After 
95 iterations, the square object is leaving the image frame 
from the far right corner. The sections of the square object 
that have left the image frame are shown reentering the image 
frame from the appropriate, opposite edges. 

The EKF was able to detect the square object motion as it 
traversed the homogeneous background. The high noise level 
caused the EKF to take considerably longer to estimate the 
object velocity, but the EKF was still able to provide an 
enhance estimate of the image frame (Figure 4.17). 

The increased noise level and addition of the checkerboard 
background precluded the EKF from detecting the square test 
object and estimating its velocity. 

The next set of simulations depicted in Figures 4.24 
through 4.31 were included to show how the EKF performed when 
the standard deviation of the noise was increased to 4.0. 
After 100 iterations, the EKF was unable to detect the square 
object or provide a reasonable estimate of the object 
velocity. 

The square object velocities of simulations #7 through 
#10 included noninteger values. The standard deviation of 
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the noise was kept at 1.0 for all four simulations. The 
response and performance of the EKF was similar to the first 
two simulations where the square object moved with integer 
velocities and the standard deviation of the noise was 1.0. 

The computer simulation results using the pyramid test 
object are depicted in Figures 40 through 71. The various 
test conditions for the pyramid object are summarized in Table 
4.2. 

The pyramid test object contains slightly less than half 
of the signal energy content in the square test object. This 
had a direct effect on the performance of the EKF. 

Simulations #1 and #2 are shown in Figures 4.48 through 
4.55. The first two simulations using the pyramid test object 
had identical test conditions as the first two simulations 
using the square test object. However, the EKF needed about 
25 more iterations before it was able to track the motion of 
the pyramid test object (Figure 4.50). 

The addition of the checkerboard background in simulation 
#2 further reduced the signal energy of the pyramid test 
object. This adversely affected the performance of the EKF 
as it was unable to detect or track the motion of the pyramid 
test object. 

In simulations #3 and #4, the standard deviation of the 
noise was raised to 2.0. Here, the increased noise level 
prevented the EKF from detecting the pyramid test object. 
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TABLE 4.2 SIMULATION TEST CONDITIONS (PYRAMID OBJECT) 



Simulation 


Background 


Velocity 


Noise level 


1 


H 


(1.0, 2.0) 


1.0 


2 


C 


(1.0, 2.0) 


1.0 


3 


H 


(1.0, 2.0) 


2.0 


4 


C 


(1.0, 2.0) 


2.0 


5 


H 


(2. 5, 1.5) 


1.0 


6 


C 


(2.5, 1.5) 


1.0 


7 


H 


(1.0, 0.5) 


1.0 


8 


C 


(1.0, 0.5) 


1.0 



H - homogeneous background, C - checkerboard background 



Simulations #5 through #8 had the pyramid test object 
moving at noninteger velocities. The standard deviation of 
the noise was kept at 1.0 for each of these simulations. The 
EKF was able to detect and closely estimate the noninteger 
velocities for simulations #5 and #7 which had the pyramid 
test object moving across the homogeneous background. Again, 
as was experienced in simulation #2, the addition of the 
checkerboard background precluded the EKF from detecting and 
tracking the pyramid test object. 
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for Simulation #1 after 40 Iterations. 
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Figure 4.10 Velocity estimates for Simulation #1. 




Iteration 



Figure 4.11 Frobenius norm results for Simulation #1. 
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Figure 4.12 Image Frame for Simulation #2 after 40 
Iterations. 




Figure 4 . 13 Extended Kalman Filter Estimate of Current 
Image Frame for Simulation #2. 






Figure 4.14 Veclocity estimates for Simulation #2. 




Figure 4.15 Frobenius norm results for Simulation #2. 
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Figure 4.18 Velocity estimates for Simulation #3. 
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Figure 4.20 
Iterations. 



Image Frame for Simulation #4 after 95 
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Figure 4.22 Velocity estimates for Simulation #4. 




Iteration 

Figure 4.23 Frobenius norm results for Simulation #4. 
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Iterations. 




Frame for Simulation #5 after 100 Iterations. 
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Figure 4.26 Velocity Estimates for Simulation #5. 




Figure 4.27 Frobenius Norm Results for Simulation #5. 
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Iterations. 




Frame for Simulation #6 after 100 Iterations. 

t 
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Figure 4.30 Velocity Estimates for Simulation #6. 




Figure 4.31 Frobenius Norm Results for Simulation #6. 
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Figure 4.32 
Iterations. 



Image Frame for Simulation #7 after 45 




Figure 4.33 Extended Kalman Filter Estimate of Image 
Frame for Simulation #7 after 45 Iterations. 
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Figure 4.34 Velocity Estimates for Simulation #7. 




Figure 4.35 Frobenius Norm Results for Simulation #7. 
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Iterations. 




Frame for Simulation #8 after 45 Iterations. 
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Figure 4.38 Velocity Estimates for Simulation #8. 




Figure 4.39 Frobenius Norm Results for Simulation #8. 
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Figure 4.40 
Iterations. 



Image Frame for Simulation #9 after 65 




Frame for Simulation #9 after 65 Iterations. 
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Figure 4.42 Velocity Estimates for Simulation #9. 




Figure 4.43 Frobenius Norm Results for Simulation #9. 
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Figure 4.44 Image Frame for Simulation #10 after 65 
Iterations. 




Frame for Simulation #10 after 65 Iterations. 
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Figure 4.46 Velocity Estimates for Simulation #10. 




Figure 4.47 Frobenius Norm Results for Simulation #10. 
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Frame for Simulation #1 after 40 Iterations. 
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Figure 4.51 Frobenius Norm Results for Simulation #1. 
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Iterations. 




Frame for Simulation #2 after 40 Iterations. 
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Figure 4.54 Velocity Estimates for Simulation #2. 




Figure 4.55 Frobenius Norm Results for Simulation #2. 
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Iterations. 




Frame for Simulation #3 after 65 Iterations. 
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Figure 4.58 Velocity Estimates for Simulation #3, 




Figure 4.59 Frobenius Norm Results for Simulation #3. 
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Figure 4.60 Image Frame for Simulation #4 after 65 
Iterations. 




Frame for Simulation #4 after 65 Iterations. 






Figure 4.62 Velocity Estimates for Simulation #4. 
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Iterations . 




Frame for Simulation #5 after 95 Iterations. 
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Figure 4.66 Velocity Estimates for Simulation #5. 




Figure 4.67 Frobenius Norm Results for Simulation #5. 



t 
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Iterations. 




Image Frame for Simulation #6. 
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Figure 4.70 Velocity Estimates for Simulation #6. 




Figure 4.71 Frobenius Norm Results for Simulation #6. 
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Iterations. 




Frame for Simulation #7 after 65 Iterations. 



65 





Figure 4.74 Velocity Estimates for Simulation #7. 




Figure 4.75 Frobenius Norm Results for Simulation #7. 
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Iterations. 




Frame for Simulation #8 after 65 Iterations. 
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Figure 4.78 Velocity Estimates for Simulation #8. 




Figure 4.79 Frobenius Norm Results for Simulation #8. 
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V. CONCLUSIONS 



The extended Kalman filter (EKF) , as applied to the moving 
object model, was evaluated under a variety of test 
conditions. These test conditions included various noise 
levels, two types of backgrounds, and integer and noninteger 
velocities. Inclusion of the Frobenius norm provided an 
analytical measure of the EKF's performance. The performance 
of the EKF was directly linked with its ability to detect and 
track the test object's motion. 

For simulations using a homogeneous background, the EKF 
was unable to detect and track the square test object when 
the standard deviation of the noise was raised to 4.0. 
Similar results occurred when the pyramid test object was used 
and the standard deviation was increased to 2.0. 

The detection of the test object by the EKF is also 
related to the signal energy content of the test object in 
relation to the energy contained in the image frame 
background. The addition of the checkerboard background 
increased the energy content of the background. As a result, 
the simulations showed that the detection threshold of the EKF 
was reduced. The EKF was precluded from detecting and 
tracking the square test object's motion across a checkerboard 
background when the standard deviation of the noise was raised 
above 2.0. Even with a standard deviation of 1.0, the EKF was 
unable to detect and track the pyramid test object. 
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These computer simulations have demonstrated that the EKF 



algorithm can successfully operate in a 
environment. These results are encouraging 
continued research. 



high-noise 
and warrant 
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APPENDIX A. TUTORIAL ON TWO-DIMENSONAL SPATIAL FREQUENCIES 



A. PROGRAM DISCUSSION 

The program SPAPICT.M begins by placing an 8x8 square 
object of unity amplitude in the center of the frame (Fig. 
A.l). The image frame is then transformed to the frequency 
domain where each spatial frequency component is defined by 
the discrete Fourier transform [ Ref. 2 ] . 



where 0 < ki < Ni-1 and 0 < k 2 < N 2 -I. 

The transformed image frame is sent through an ideal 
lowpass filter (LPF) , i.e., the high frequency components are 
truncated and set equal to zero. The filtered image frame is 
then transformed back to the spatial domain where the sharp 
edges are lost due to the high frequency components being 
filtered out (Fig. A. 2). 

The image frame is transformed into a matrix that exhibits 
conjugate symmetry. Conjugate symmetry results when a real 
image is transformed into the frequency domain. The number 
of computations for each spatial frequency and its complex 
conjugate can be reduced using Euler's identity, 




(A.l) 



e^“ + e"^^^ = 2 cos(w) , 



(A. 2) 
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the Spatial Domain. 
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where w = kini+k 2 n 2 for the two-dimensional case. The 
contribution of each spatial frequency and its complex 
congugate was computed using 



F(k,,k2) + F*(k,,k2) = 

2 I F(k, ,k2> I cos[^(n,k,+n2k2) +iF(k,,k2)] , (A. 3) 
where N = 32. 

As an example, the spatial frequency F(-l,-l) is the 
complex conjugate of the spatial frequency F(l,l). In the 
reconstruction of the filtered frame in the spatial frequency 
domain, the contribution of a specific spatial frequency and 
its complex conjugate were calculated using Eqn. A. 3 and then 
summed with the previous frame. 

Figures A.3-A.8 show the successive summation of each 
spatial frequency and its corresponding complex conjugate in 
the reconstruction of the filtered frame in the spatial 
frequency domain. Each figure depicting a specific spatial 
frequency has been scaled and phase shifted based on the 
Fourier coefficients that were computed using the DFT. 

Figure A. 3 depicts the summation of the dc component and 
F(0,1) and its complex conjugate. Figure A. 4 depicts the 
contribution of spatial frequencies F(0,2) and F*(0,2) which 
are then summed with the frame shown in Fig. A. 3 to produce 
Fig. A. 5. 
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Figure A. 3 Summation of the do component, F(0,1), and 
F*(0,1) . 
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Figure A. 5 Summation of Previous Spectral Frequencies, 
F(0,2) , and F*(0,2) . 
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This process is continued with the contributions from the 
other spatial frequencies being computed and summed with the 
frame shown in Fig. A. 5. Figure A. 6 shows the spatial 
frequency F(2,l) which has two cycles in the n^-direction and 
one cycle in the n 2 ~direction. Figure A. 7 shows the summation 
of all the contributions up to this point. The final frame 
shown in Fig. 8 is the same as the frame depicted in Fig. A. 2. 
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B. PROGRAM LISTING 



SPAPIC.M is a stand-alone MATLAB program that 
reconstructs a filtered image in the frequency domain. 

icicicicrkicicicicic'k'k’k’kic'kicicrk'kicrk'kicicic'k'k’kicic’k'k’k'k'kicrk'k’k’kicic’k’krk'k’k’kic’kic’ik'k'kitieitic 

% VARIABLE DEFINITION 
% 

% pict 

% 

% img 
% freq 
% 

% 

% 

% fftpic 
% 

% freqpict 
% 

% pictifft 
% 

% 

^****** ******************* ********************************** 

% CREATE THE PICTURE 

img = ones (8) ; 
pictl = zeros(32); 
pictl([14:21] , [14:21]) = img; 

[ nl n2 ] = size (pictl) ; 

% CREATE FREQUENCY VECTOR 

n = 0; 

for i = 0:3, 
if i == 0, 

for j = 0:3, 

n = n + 1; 

freq(n,:) = [ (2*pi*i)/nl (2*pi*j)/n2 ]; 

end 

else 

for j = -3:3, 

n = n + 1; 

freq(n,:) = [ (2*pi*i)/nl (2*pi*j)/n2 ]; 

end 

end 

end 



: Frame with square image as it would appear in 
spatial domain 

: 8x8 square image with unity amplitude 

: Frequency vector listing the pertinent 
spatial frequencies needed to 
reconstruct the filtered image in the 
spatial frequency domain 

; Frame after transformation to the spatial 
frequency domain 

: Frame after lowpass filtering in spatial 
frequency domain 

: Lowpass filtered Frame after transformation 
back to spatial domain 
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% PLOT THE FRAME IN THE SPATIAL DOMAIN 



mesh (pictl) 

title ('Frame — Spatial Domain') 
met a spaplot 

% TAKE FOURIER TRANSFORM OF FRAME 
fftpictl = fftshift(fft2 (pictl) ) ; 

% PASS THE TRANSFORMED FRAME THROUGH A LOWPASS FILTER 

trnpict - f ftpictl ( [14 : 20] , [ 14 : 20] ) ; 
trnpictl = trnpict(:); 

% TAKE INVERSE FOURIER TRANSFORM OF FILTERED FRAME 

freqpict = zeros(32); 

freqpict( [14:20] , [14:20])= trnpict; 

pictifft = if ft2 (fftshift ( freqpict) ) ; 

% PLOT THE FILTERED FRAME IN THE SPATIAL DOMAIN 

mesh (pictifft) 

title (' Lowpass Filtered frame in the spatial domain') 
meta spaplot 
pause (5) 

% DISPLAY AND SAVE THE MAGNITUDE AND PHASE OF THE % 
TRANSFORMED PICTURE 

diary spadiary.mat 

disp (' Lowpass Filtered Frame — Magnitude and Phase') 
zl = [ abs (trnpictl' ) ; angle (trnpictl ' ) ] 
diary off 

% CONSTRUCT THE FRAME IN THE FREQUENCY DOMAIN USING 
% FREQUENCIES THAT WERE NOT TRUNCATED 

Zl = zl(: ,1:25) ; 
zl = flipud(zl'); 
pict = zeros(32); 

for n = 1:25, 

u = freq(n, 1) ; 

V = freq(n, 2) ; 

mag = z 1 ( n , 1 ) ; 

phase = zl (n, 2) ; 

X =0:31; 

y =0:31; 

[ XX, yy ] = meshdom(x,y) ; 

z = 2 * mag * cos(u*xx + v*yy + phase) ; 
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pict = pict + z; 
subplot (211) , mesh(z) 

title ( 'Constructing truncated image in frequency 
'domain' ) 

xlabel ( [ ' f ( ' ,num2str (u) , ' , ' , num2str (v) ,')']) 
subplot (212) , mesh (pict) 

xlabel ([ 'Adding f ( ' ,num2str (u) , ' , ' ,num2str (v) , ' ) and'... 
' f*(', num2str (u) , ' , ' , num2str (v) , ' ) to the '... 
'picture' ] ) ; 
meta spaplot 
pause (3) 
clg 

end 

clg, subplot(lll) 

% GET A PLOT OF THE FINAL IMAGE CONSTRUCTED IN THE FREQUENCY 
% DOMAIN 

mesh (pict) 

title( 'Final image constructed in the frequency domain') 
meta spaplot 
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APPENDIX B. PPHSE.M 



A. PROGRAM DISCUSSION 

PPHSE.M is a stand-alone program written for use with 
MATLAB to demonstrate the spatial frequency phase shift that 
results from moving an object across an image frame in the 
spatial domain. 

PPHSE.M sequentially moves an 8x8 unity square across a 
32x32 image frame in the spatial domain for 50 iterations. 
Two image frames are transformed to the spatial frequency 
domain. The first image frame consists of a homogeneous 
background with no added noise. The second image frame has 
zero mean, Gaussian white noise. 

The real and imaginary Fourier coefficients and the phase 
in degrees are recorded and plotted for both the no-noise and 
noise-corrupted image frames for the two spatial frequencies 
F(1,0) and F(0,1). Only the phase is worth noting since the 
magnitude of the transformed image frame remains constant. 
In the spatial frequency domain, the effect of moving the 
square object across the image frame is a phase shift in each 
of the spatial frequency components. 

Figure B.l shows the phase shift in degrees for the spatial 
frequency F(1,0). The effect of moving the square object 
across the image frame 1 pixel per iteration is a 11.25° phase 
shift per iteration. The dashed line represents the noise- 
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corrupted image frame. The phase shift is no longer perfectly 
linear, but it still tracks the no-noise phase shift rather 
closely. 

Figure B.2 shows the phase shift in degrees for the 
spatial frequency F(0,1) . The effect of simultaneously moving 
the square object 2 pixels per iteration vertically is a 22.5° 
phase shift per iteration. The phase shift for the noise 
corrupted image frame is represented by the dashed line and 
the nonlinear phase shift again generally follows the no-noise 
phase shift. 

It should be noted that the phase shift due to the object 
motion is a function of the spatial frequency. In this 
example, F(1,0) and F(0,1) were chosen because they represent 
the spatial frequencies that purely describe the horizontal 
and vertical motion of the moving object. 

Figures B.3-B.6 sequentially show the magnitude of the 
real and imaginary Fourier coefficients for the spatial 
frequencies F(1,0) and F(0,1). The dashed line again 
represents the coefficients of the noise-corrupted image 
frame. 

The real and imaginary Fourier coefficients for each of 
the spatial frequencies were two of the three the states that 
were estimated by the extended Kalman filter described in this 
thesis . 
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Figure B.l Phase Shift in Degrees for Spatial Frequency 
F(1,0) . 




Figure B.2 Phase Shift in Degrees for Spatial Frequency 
F(0,1) . 
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Figure B.3 Real Fourier Coefficients for Spatial 
Frequency F ( 1 , 0 ) . 




Figure B.4 Imaginary Fourier Coefficients for Spatial 
Frequency F(1,0). 



/ 
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Figure B.5 Real Fourier Coefficients for Spatial 
Frequency F ( 0 , 1 ) . 




Figure B.6 Imaginary Fourier Coefficients for Spatial 
Frequency F(0,1). 



t 
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B. PROGRAM LISTING 



% PPHSE.M is a stand-alone program that demonstrates 
% the spatial frequency phase shift that results when 
% a moving object traverses an image frame in the 
% spatial domain 

clear % clear out all variables in memory 

clear functions 



shpe = 
vrow = 
vclm = 
imax = 
jmax = 
its = 
sigl = 



['Square']; % Define the shape of the moving image 
2.0; % Vertical Velocity — # of pixels/ its 

1.0; % Horizontal Velocity — # of pixels/its 

3; % Size of the spatial frequency truncation 

3; 

50; % Number of iterations 

1.0; % Standard deviation of the noise 



shape = ones (8); 
pict = zeros(32); 



% Define the shape — square of % 
unity amplitude 

% Create digital array to place % 
image in 

rand( 'normal' ) ; % Establish normal distribution 

rand( 'seed' , 2380849) % Define a standard random seed 

[ nl n2 ] = size (shape) ;% Define the size of the image 



% CREATE STORAGE VECTORS FOR CALULATED DATA 



phase0=zeros (2 , its) ; % Phase shift in degrees — no noise 
phasel=zeros (2 , its) ; % Phase shift in degrees — 1.0 Std dev 

rslto = zeros (2 , its) ; % Fourier coef. — no noise 
rsltl = zeros (2 , its) ; % Fourier coef. — Std dev. = 1.0 

% MOVE THE IMAGE ACROSS THE DIGITAL ARRAY FOR THE NUMBER OF 
% ITERATIONS SPECIFIED. 



for i = 0: (its-1) , 

% DEFINE THE PARAMETERS FOR IMAGE LOCATION 

11 = rem(round(vrow*i+l) , 32)+l; 

12 = rem ( round (vrow*i+8 ), 32) +1; 

13 = rem ( round (vclm*i+l ), 32) +1; 

14 = rem (round (vclm*i+8) , 32) +1 ; 
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% PUT IMAGE INTO THE 32X32 ARRAY "pict" 

pict = zeros (32) ; 

if i2 < il, 
if i4 < i3, 



pict(il:32, 13:32) = 
pict ( 1: i2 , 1 : i4) = 

pict (il: 32 , 1 : i4) = 

pict(l:i2,i3:32) = 

else 

pict (il: 32 , i3 : i4) = 
pict (1: i2 , i3 : i4) = 

end 

else 

if i4 < i3, 

pict(il:i2,i3:32) = 
pict(il: i2, l:i4) = 

else 

pict (il: i2 , i3 : i4) = 

end 

end 



shape (1: nl-i2, I:n2-i4) ; 
shape(nl+l-i2:nl,n2+l-i4 :n2) ; 
shape ( 1 : nl-i2 , n2+l-i4 : n2 ) ; 
shape (nl+l-i2 :nl, I:n2-i4) ; 

shape (1 : nl-i2, l:n2) ; 
shape(nl+l-i2:8, l:n2) ; 



shape ( 1 :nl, I:n2-i4) ; 
shape ( 1 : nl , n2+l-i4 : n2 ) ; 

shape ( 1 : nl , 1 : n2 ) ; 



% CREATE THE FRAME CORRUPTED BY NOISE 



pictl = pict + (sigl*rand (32) ) ; 

% TRANSFORM THE FRAMES INTO THE SPATIAL FREQUENCY DOMAIN 



fftpictO = fftshift(fft2 (pict) ) ; 
fftpictl = fftshift(fft2 (pictl) ) ; 



% OBTAIN THE PHASE IN DEGREES FROM SPATIAL FREQUENCIES 
% F(0,1) & F(1,0) 



phaseO ( 1 , i+1) = 57 . 2958*angle (fftpictO (16, 17) ) ; 
phasel (1, i+1) = 57.2958*angle(fftpictl(16, 17) ) ; 
phaseO (2 , i+1) = 57.2958*angle(fftpict0 (17, 18) ) ; 
phasel(2,i+l) = 57.2958*angle(fftpictl(17, 18) ) ; 



% OBTAIN THE FOURIER COEFFICIENTS FOR SPATIAL FREQUENCIES 
% F(0,1) & F(1,0) 



rslto ( 1 , i+1) 
rsltl (1, i+1) 
rslto (2 , i+1) 
rsltl (2 , i+1) 



fftpictO (16,17) ; 
fftpictl (16, 17) ; 
fftpictO (17, 18) ; 
fftpictl (17, 18) ; 



end 
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% PLOT THE PHASE SHIFT FOR SPATIAL FREQ F(0,1) IN DEGREES 
n = l:its; 

plot (n,phaseO (1, : ) , 'w' ,n,phasel(l, : ) , ' — w' ) 
title ('Phase shift for image velocity of (1,2)') 
xlabel ( 'Phase shift for F(0,1) in degrees') 
ylabel ( ' Degrees ' ) 
meta pphse 
pause 

% PLOT THE PHASE SHIFT FOR SPATIAL FREQUENCY F(1,0) IN 
% DEGREES 

plot (n,phaseO (2 , : ) , 'w' , n,phasel (2 , ; ) , ' — w' ) 
title ('Phase shift for image velocity of (1,2)') 
xlabel ( 'Phase shift for F(1,0) in degrees') 
ylabel ( ' Degrees ' ) 
meta pphse 
pause 

% PLOT THE REAL FOURIER COEFF FOR SPATIAL FREQ F(0,1) 

plot (n, real (rsltO (1 , : ) ) , 'w' , n, real (rsltl (1, : ) ) , ' — w' ) 
title('Real Fourier coefficients for F(0,1)') 
xlabel (' Iteration' ) , ylabel ( 'Magnitude' ) 
meta pphse 
pause 

% PLOT THE IMAGINARY FOURIER COEFF FOR SPATIAL FREQ F(0,1) 

plot (n, imag (rsltO (l,:)),'w',n, imag (rsltl (1 , : ) ) , ' — w' ) 
title (' Imaginary Fourier coefficients for F(0,1)') 
xlabel ( ' Iteration' ) , ylabel ( 'Magnitude' ) 
meta pphse 
pause 

% PLOT THE REAL FOURIER COEFF FOR SPATIAL FREQ F(1,0) 

plot (n, real (rsltO (2 , : ) ) , 'w' , n, real (rsltl (2 , : ) ) , ' — w' ) 
title('Real Fourier coefficients for F(1,0)') 
xlabel ( ' Iteration' ) , ylabel ( 'Magnitude' ) 
meta pphse 
pause 

% PLOT THE IMAGINARY FOURIER COEFF FOR SPATIAL FREQ F(1,0) 

plot (n, imag (rsltO (2 , : ) ) , 'w' , n, imag (rsltl (2 , : ) ) , ' — w' ) 
title (' Imaginary Fourier coefficients for F(1,0)') 
xlabel ( 'Iteration' ) , ylabel ( 'Magnitude' ) 
meta pphse 
pause 
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APPENDIX C. PMEKF.M 



PROGRAM LISTING 

PMEKF.M is a stand-alone MATLAB program that implements the 

modified extended Kalman filter (MEKF) for moving objects. 

********************************************************** 
VARIABLE DEFINITIONS 



shpe 

vrow 

vclm 

imax, jmax 

its 

sig 

chkrbd 

a 

runnum 

pictchkr 

freq 

Ifreq 

normpic 

normpict 



xh 

vh 

y 

pk 

pkk 

q 

r 

pict 

pictn 

randpict 

********** 



: string defining shape of moving object 
: vertical velocity in # of pixels per iteration 
: horiz. velocity in # of pixels per iteraiton 
: rectangular dimensions of spatial truncation 
window; located at center of transformed frame 
: number of iterations 

: standard deviation of the added noise 
: flag indicating user's intentions of whether 
or not to include checkerboard background 
: amplitude of individual checkers 
: for successive runs; number indicating 
specific run 

: image frame containing checkerboard bkgrd; 

moving image included later 
: vector listing two-dimensional frequencies to 
be analyzed 

: length of the freq vector 
: Frobenius norm of error matrix; subtract 
filtered image frame from original, no-noise 
image frame 

: Frobenius norm of error matrix; subtract 

filtered image frame from truncated, no-noise 
image frame 

; state estimate vector 
: velocity estimate vector 
: vector containing Fourier coefficients 
: storage array for **pkk“ covariance matrices 
: estimation error covariance matrix; spatial 
frequency specific 
: plant noise covariance matrix 
: measurement noise covariance matrix 
: no-noise image frame containing moving object 
: no-noise image frame containing moving object 
and spatial frequency truncated 
: image frame containing moving object and added 
zero-mean, white Gaussian noise 
*********************************************** 
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clear 

clear functions 



[ 'Rectangle' ] ; 

2 . 0 ; 

1 . 0 ; 

3; 



shpe = 
vrow = 
vclm = 
imax = 
jmax = 
its = 
sig = 
chkrbd 
a = 
runnum 



3; 

25; 

1 . 00 ; 

= 'n'; 

0.25; 

= i; 



% ( 'y' 



% DEFINE THE SHAPE 



yes; 'n' 



no ) 



xl = 1. 


o 

o 


X2 


= 1.00; x3 = 


= 1, 


.00; 


; x4 


= 1.00; 


shape = 


[[ 


xl 


xl 


xl 


xl 


xl 


xl 


xl 


xl 


] 




[ 


xl 


x2 


x2 


x2 


x2 


x2 


x2 


xl 


] 




[ 


xl 


x2 


x3 


x3 


x3 


x3 


x2 


xl 


] 




[ 


xl 


x2 


X3 


X4 


X4 


X3 


x2 


xl 


] 




[ 


xl 


x2 


x3 


x4 


x4 


x3 


x2 


xl 


] 




[ 


xl 


x2 


X3 


X3 


X3 


X3 


x2 


xl 


] 




[ 


xl 


x2 


x2 


x2 


x2 


x2 


x2 


xl 


] 




[ 


xl 


xl 


xl 


xl 


xl 


xl 


xl 


xl 


]]; 



% CREATE A CHECKERBOARD BACKGROUND IF REQUIRED 

i f ( chkrbd == ' y ' ) , 

pictchkr = zeros (32) ; 
chkr = a * ones ( 8 ) ; 

for i = 1:4, 

if i == 1 I i ==3, 

pictchkr ([ i*8-7 : i*8 ], [1 : 8 ] ) = chkr; 
pictchkr ( [i*8-7 : i*8] , [17:24]) = chkr; 
else 

pictchkr ( [i*8-7 : i*8] , [9 : 16] ) = chkr; 
pictchkr ([ i*8-7 : i*8] , [25:32]) = chkr; 

end 

end 

end 
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% CREATE FREQUENCY VECTOR 
n = 0; 



for i = 0:imax, 
if i == 


0, 






for 


j = 0:jmax, 


else 


end 

for 


n = n+1; 
freq(n,:) = 

j = -jmaxtjmax 


end 

end 


end 


n = n+1; 
freq(n,:) ^ 



% CREATE AND DEFINE NECESSARY VARIABLES 

CO = pi/16; 
j j = sqrt(-l) ; 

Ifreq = length (freq) ; 

rand ( 'normal ') ; % Establish normal distribution 

rand( 'seed' , 2380849) % Define a standard random seed 

[ nl n2 ] = size (shape) ;% Dimensions of the image frame 
y = zeros (Ifreq, its) ; 
xh = zeros (3 , Ifreq) ; 
vh = zeros (2 , its) ; 
q = [ 0 0 0 

0 0 0 
0 0.4]; 

r = sig^2*512*eye(2) ; 



normpic = zeros (1 , its) ; 
normpict= zeros (1 , its) ; 

% GENERATE THE INITIAL ESTIMATION ERROR COVARIANCE MATRIX. 

% EACH INDIVIDUAL FREQ COVARIANCE MATRIX IS STACKED IN "pk" 

pk = zeros (9 , Ifreq) ; 
for k = 1: Ifreq, 

pkk = [ 1024 0 0 

0 1024 0 

0 0 4*freq(k, : ) *freq(k, ; ) ' ] ; 

pk( : ,k)= pkk( : ) ; 

end 
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% MODIFIED EXTENDED KALMAN FILTER 



% OBJECT IS MOVED ACROSS THE IMAGE FRAME FOR THE NUMBER OF 
% ITERATIONS SPECIFIED 

for i = 0:(its-l), 

% REINITIALIZE THE TYPE OF IMAGE FRAME USED IN THE 
% SIMULATION 

if ( chkrbd == ' y ' ) , 
pict = pictchkr; 

else 

pict = zeros (32) ; 

end 



% DEFINE THE PARAMETERS FOR IMAGE LOCATION 



11 = rem (round (vrow*i+l) ,32) +1; 

12 = rem(round(vrow*i+8) ,32)+l; 

13 = rem(round(vclin*i+l) ,32)+l; 

14 = rem (round (vclm* i+8 ) ,32) +1; 



% PUT IMAGE INTO THE 32X32 ARRAY "pict” 



if i2 < il, 
if i4 < i3, 

pict(il:32,i3:32) 
pict(l:i2,l:i4) 
pict (il : 32 , 1: i4 ) 
pict(l:i2,i3:32) 
else 

pict(il:32,i3:i4) 
pict (1 ; i2 , i3 : i4 ) 
end 

else 

if i4 < i3, 

pict(il:i2,i3;32) 
pict (il : i2 , 1; i4) 
else 

pict(il:i2,i3:i4) 

end 

end 



shape (1 :nl-i2, I:n2-i4) ; 
shape (nl+l-i2 : nl, n2+l-i4 : n2 ) ; 
shape ( 1 : nl-i2 , n2+l-i4 : n2 ) ; 
shape (nl+l-i2 : nl, 1 : n2-i4) ; 

shape ( 1 : nl-i2 , 1 : n2 ) ; 
shape ( nl+l-i2 : 8 , 1 : n2 ) ; 



shape (l;nl, I;n2-i4) ; 
shape ( 1 : nl , n2+l-i4 : n2 ) ; 

shape ( 1 : nl , 1 : n2 ) ; 



randpict 

fpic 

fpic 

fband 

fvec 

y(: ,i+l) 



pict + (sig*rand(32) ) ; 
fft2 (randpict) ; 
f ftshift (fpic) ; 

fpic( (17-imax) : (17+imax) , (17-jmax) : (17+jmax) ) ; 
fband ( : ) ; 

fvec( ( ( length ( fvec) +1) /2) ; length(fvec) ) ; 
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% TRUNCATE THE SPATIAL FREQS OF UNCORRUPTED PICTURE 

pictn = fft2(pict); 

pictn = fftshift (pictn) ; 

fbandl =pictn( (17-imax) : (17+imax) , (17-jmax) : (17+jmax) ) ; 
pictn = zeros(32); 

pictn ( (17-imax) : (17+iinax) , (17-jmax) : (17+jmax) ) = fbandl; 
pictn = fftshift (pictn) ; 

pictn = if ft2 (pictn) ; 

% ITERATE THRU THE "parallel" BANK OF EXT. KALMAN FILTERS 

for k = Itlfreq, 

dth = co*xh(3,k); 
cdth = cos (dth) ; 
sdth = sin (dth) ; 
dh = [ cdth sdth 

-sdth cdth ] ; 

dph = [ (-xh(l,k) *sdth+xh(2,k) *cdth) *co 

(-xh(l,k) *cdth-xh(2,k) *sdth) *co ] 

ah = [ dh dph 

0 0 1]; 

pkk( : ) = pk( : ,k) ; 
pkplk = ah*pkk*ah' +q; 
g = pkplk( ; , 1: 2)/ (pkplk(l; 2 , 1: 2) +r) ; 
xh(l:2,k) = dh*xh (1 : 2 , k) ; 
yh = xh(l;2,k) ; 
ym = [ real (y (k, i+1) ) 

imag(y(k,i+l) ) ]; 
xh(:,k) = xh ( : , k) +g* (ym-yh) ; 
pkk = pkplk - g* [1 0 0; 0 1 0]*pkplk; 
pk( : ,k) = pkk( : ) ; 

end 

THESE THREE LINES GENERATE THE WEIGHTED LEAST SQUARES 
ESTIMATE OF THE VELOCITY VECTOR. COMMENTED OUT WHEN 
NOT USED. 

sigin = diag(ones(l, (lfreq-1) ) ./pk(9 , 2 : Ifreq) ) ; 

fre = freq(2 : Ifreq, :) ; 

vh(;,i+l) = (fre^*sigin*fre)\fre'*sigin*xh(3,2;lfreq) ' ; 

% THIS LINE GENERATES THE LEAST SQUARES ESTIMATE OF THE 
% VELOCITY VECTOR (unweighted) . COMMENTED OUT WHEN NOT 
% USED. 

% vh(;,i+l) = f req(2 ; Ifreq, : ) \xh (3 , 2 ; Ifreq) ' ; 



94 



% THIS LINE CHANGES THE ESTIMATES OF FREQUENCY FOR THE 
% INDIVIDUAL FILTERS LEAST SQUARES ESTIMATE. COMMENTED OUT 
% WHEN NOT USED. 

xh(3,:) = [ freq*vh( ; , i+1) ]'; 



% RECONSTRUCT THE IMAGE BACK INTO THE TIME DOMAIN. 

fpic = zeros (32) ; 

xhc = zeros (1, (2*lfreq-l) ) ; 

xhc (Ifreq; (lfreq*2-l) ) = xh(l, : )+j j*xh(2, : ) ; 
xhc^freq;-!: 1) = xh(l, : ) -j j*xh(2 , : ) ; 
xhc(lfreq) = xh(l,l); 
f band ( : ) = xhc ; 

fpic( (17-iraax) : (17+imax) , (17-jmax) : (17+jmax) ) = fband; 
fpic = fftshift (fpic) ; 
pic = ifft2(fpic); 

% CALCULATE THE NORM OF A MATRIX 

% FROBENIUS NORM : UNCORRUPTED IMAGE 

normpic(i+l) = norm(pict - pic, 'fro')/ 32; 

% FROBENIUS NORM : UNCORRUPTED, SPATIALLY -TRUNCATED IMAGE 

normpict(i+l)= norm(pictn - pic, 'fro')/ 32; 

end 

% DEFINE A TEXT LABEL FOR FINAL PLOTS 

name =['Std deviation = ' ,num2str (sig) , ' — '... 

'Iterations = ' , num2str (its) , ' — '... 

'vel ( ' , num2str (vrow) , ' , ' ,num2str (vclm) ,')']; 

% DISPLAY THE FINAL IMAGE 

pname = ['meta pmekf_' , num2str (runnum) ]; 

clc 

mesh (pic) 

title ([ shpe ' — Filtered Image — Time domain']) 
xlabel (name) 
eval (pname) 
pause 
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mesh (pict) 

title ([ shpe ' — Uncorrupted image']) 
xlabel (name) 
eval (pname) 
pause 

mesh (pictn) 

title ([ shpe ' — Uncorrupted image (spatially truncated)']) 
xlabel (name) 
eval (pname) 
pause 

% PLOT FROBENIUS NORMS AND VELOCITY ESTIMATES 
n = l:its; 

plot ( n , normpic , ' - ' , n , normpict , ' — ' ) 
title ([ shpe ' — Frobenius Norm']) 
eval (pname) 
pause 

plot(n,vh(l, :),'-' ,n,vh(2, :), ' — ') 
title ([ shpe ' — Velocity Estimates']) 
xlabel (name) 
eval (pname) 
pause 

% SAVE THE CALCULATED DATA 
m = [ n' vh' normpic' normpict']; 

ekfdata=[ 'save pmekf_' , num2str (runnum) , ' .mat m /ascii']; 
eval (ekfdata) ; 
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